    #include <iostream>
    #include <cmath>
    #include <cstdlib>
    #include <cuda_runtime.h>
    
    // Macro for error checking of CUDA API calls
    #define CHECK_CUDA_ERROR(call) { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            std::cerr << "CUDA error at " << __FILE__ << ":" << __LINE__ << ": " \
                      << cudaGetErrorString(err) << std::endl; \
            exit(EXIT_FAILURE); \
        } \
    }
    
    // Kernel function to compute sliding window k-quantile for each window position
    __global__ void sliding_window_k_quantile_kernel(
        const double* input, 
        double* output, 
        const int rows, 
        const int cols, 
        const int window_width, 
        const int result_cols, 
        const double k
    ) {
        // Calculate global window index: 
        //   blockIdx.x - row index
        //   blockIdx.y - block index along window dimension
        //   threadIdx.x - thread index within block
        const int window_index = blockIdx.y * blockDim.x + threadIdx.x;
        if (window_index >= result_cols) return;  // Check if window index is within bounds
    
        const int row = blockIdx.x;  // Current row index
    
        // Each thread handles one window: load window elements into local array
        double window[64];  // Static array (max window_width=64)
        const int start_col = window_index;
        for (int i = 0; i < window_width; ++i) {
            window[i] = input[row * cols + start_col + i];
        }
    
        // Insertion sort for the window (efficient for small arrays)
        for (int i = 1; i < window_width; ++i) {
            double key = window[i];
            int j = i - 1;
            // Shift elements greater than key to the right
            while (j >= 0 && window[j] > key) {
                window[j + 1] = window[j];
                --j;
            }
            window[j + 1] = key;
        }
    
        // Compute quantile using linear interpolation (matching numpy.quantile)
        const double index = k * (window_width - 1);
        const int low_index = static_cast<int>(floor(index));
        const int high_index = static_cast<int>(ceil(index));
        const double fraction = index - low_index;
    
        double quantile;
        if (low_index == high_index) {
            quantile = window[low_index];
        } else {
            quantile = window[low_index] + fraction * (window[high_index] - window[low_index]);
        }
    
        // Write result to output matrix (row-major order)
        output[row * result_cols + window_index] = quantile;
    }
    
    int main() {
        // Example configuration (matching the Python example)
        const int rows = 1024;
        const int cols = 1024;
        const int window_width = 16;
        const double k = 0.5;
    
        // Validate window_width
        if (window_width > 64) {
            std::cerr << "Error: Window width exceeds maximum supported size (64)." << std::endl;
            return EXIT_FAILURE;
        }
    
        // Allocate host memory for input and output
        double* host_matrix = new double[rows * cols];
        const int result_cols = cols - window_width + 1;
        double* host_result = new double[rows * result_cols];
    
        // Initialize input matrix with random values (seed=42)
        std::srand(42);
        for (int i = 0; i < rows * cols; ++i) {
            host_matrix[i] = static_cast<double>(std::rand()) / RAND_MAX;
        }
    
        // Allocate device memory
        double* d_input, * d_output;
        CHECK_CUDA_ERROR(cudaMalloc(&d_input, rows * cols * sizeof(double)));
        CHECK_CUDA_ERROR(cudaMalloc(&d_output, rows * result_cols * sizeof(double)));
    
        // Copy input data to device
        CHECK_CUDA_ERROR(cudaMemcpy(d_input, host_matrix, rows * cols * sizeof(double), cudaMemcpyHostToDevice));
    
        // Configure kernel launch parameters
        const int block_size = 256;  // Threads per block
        const int blocks_per_row = (result_cols + block_size - 1) / block_size;  // Blocks per row dimension
        dim3 grid_dim(rows, blocks_per_row);  // Grid: rows x blocks_per_row
        dim3 block_dim(block_size);           // Block: block_size threads
    
        // Launch kernel
        sliding_window_k_quantile_kernel<<<grid_dim, block_dim>>>(
            d_input, d_output, rows, cols, window_width, result_cols, k
        );
        CHECK_CUDA_ERROR(cudaGetLastError());  // Check for kernel launch errors
        CHECK_CUDA_ERROR(cudaDeviceSynchronize());  // Wait for kernel completion
    
        // Copy results back to host
        CHECK_CUDA_ERROR(cudaMemcpy(host_result, d_output, rows * result_cols * sizeof(double), cudaMemcpyDeviceToHost));
    
        // Output sample information (avoid full print for large matrices)
        std::cout << "Input shape: " << rows << "x" << cols << std::endl;
        std::cout << "Output shape: " << rows << "x" << result_cols << std::endl;
        std::cout << "Sample output element [0][0]: " << host_result[0] << std::endl;
    
        // Cleanup
        delete[] host_matrix;
        delete[] host_result;
        cudaFree(d_input);
        cudaFree(d_output);
    
        return 0;
    }